############## 
############## 
############## Models
remove(list = ls())

base::library(conflicted)
base::library(tidyverse)
conflict_prefer("filter","dplyr")
base::library(ggplot2)
base::library(dplyr)
base::library(here)
conflict_prefer("here", "here")
base::library(systemfit)
#base::library(sf)
#base::library(sp)
#base::library(spdep)
#base::library(spatialreg)
#base::library(rgdal)
#base::library(utils)
#base::library(geodist)
base::library(stargazer)
#base::library(texreg)
base::library(DataCombine)
#base::library(reshape2)

load(here("Data", "Model_Data.RData"))
glimpse(final_data)

summary(final_data$death_rate)
table(final_data$year)

final_data <- final_data %>% 
  mutate(time_trend = ifelse(year == 2004, 1, NA),
         time_trend = ifelse(year == 2008, 2, time_trend),
         time_trend = ifelse(year == 2012, 3, time_trend),
         time_trend = ifelse(year == 2016, 4, time_trend),
         time_trend = ifelse(year == 2020, 5, time_trend)) 

glimpse(final_data)

eq1_d <- dem_rep ~ death_rate + unemp_rate + 
  log(defl_pcincome) + 
  share_black + share_hisp + share_asian + share_young + share_elder + 
  log1p(share_college) + log1p(share_manuf) +
  log(tot_pop) + as.factor(rural_urban) + state_partiesdiff +
  state + rep_inc + incumbent_cand + lag_dem_rep + time_trend
eq2_d <- abs_rep ~ death_rate + unemp_rate +  
  log(defl_pcincome) + 
  share_black + share_hisp + share_asian + share_young + share_elder + 
  log1p(share_college) + log1p(share_manuf) +
  log(tot_pop) + as.factor(rural_urban) + state_partiesdiff +
  state + rep_inc + incumbent_cand + lag_abs_rep + time_trend
Syst_death <- list(eq1_d, eq2_d)
model_d <- systemfit(Syst_death, method = "SUR", data = final_data)
summary(model_d)


summary(final_data$death_rate)
test_data <- final_data %>% filter(death_rate <= 5) %>% 
  select(change_death, death_rate, lag_death_rate)
summary(test_data$change_death)
summary(final_data$death_rate)
hist(final_data$death_rate)
3*sd(final_data$death_rate, na.rm = T)
#10.68225

cor(final_data$rep_share2, final_data$death_rate, use = "complete.obs")
cor(final_data$dem_share2, final_data$death_rate, use = "complete.obs")
cor(final_data$eff_abstention, final_data$death_rate, use = "complete.obs")

################################################################################
################################################################################
################################################################################
###################################Bootstrap#################################### 
################################################################################
################################################################################
################################################################################

sum(is.na(final_data$death_rate))
final_data <- final_data %>% filter(!is.na(death_rate)) 
plot(density(final_data$death_rate))

final_data2 <- final_data %>% 
  mutate(death_rate = death_rate+sd(final_data$death_rate))
glimpse(final_data2)


rep <- length(final_data$dem_rep)
indices_rep <- 1:rep
#Republican Party, Democratic Party, Abstention
results <- as.data.frame(matrix(data = NA, nrow = 3000, ncol = 3))

sample(1:1000000, 1)
set.seed(822535)

#Moving from 5 to 85 deaths

for(i in 1:3000) {
  
  # Resampling w/ replacement
  tryCatch({sample_indices <- sample(indices_rep, rep, replace = T)
  bootstrap_data <- final_data[sample_indices, ]
  
  
  #SUR Model: Republican as the baseline category
  model <- systemfit(Syst_death, method = "SUR", data = bootstrap_data)
  
  yhat <- predict(object = model, newdata = bootstrap_data) %>% 
    #filter(!is.na(eq1.pred)) %>% 
    rename(dem = eq1.pred,
           abst = eq2.pred)
  
  yhat <- yhat %>% 
    mutate(rephat = 1/(1+exp(dem)+
                         exp(abst))) %>% 
    mutate(demhat = exp(dem)/(1+exp(dem)+
                                exp(abst))) %>% 
    mutate(absthat = exp(abst)/(1+exp(dem)+
                                  exp(abst))) 
  
  yhat <- yhat %>% 
    select(rephat:absthat)
  
  ############################################################################## 
  ############################################################################## 
  ##############################################################################
  #80
  #75
  #SUR Model 1SD shock: Republican as the baseline category
  bootstrap_data2 <- bootstrap_data %>% 
    mutate(death_rate = death_rate+60)
  
  yhat2 <- predict(object = model, newdata = bootstrap_data2) %>% 
    #filter(!is.na(eq1.pred)) %>% 
    rename(dem = eq1.pred,
           abst = eq2.pred)
  
  yhat2 <- yhat2 %>% 
    mutate(rephat2 = 1/(1+exp(dem)+
                          exp(abst))) %>% 
    mutate(demhat2 = exp(dem)/(1+exp(dem)+
                                 exp(abst))) %>% 
    mutate(absthat2 = exp(abst)/(1+exp(dem)+
                                   exp(abst))) 
  
  yhat2 <- yhat2 %>% 
    select(rephat2:absthat2)
  
  ############################################################################## 
  ############################################################################## 
  ##############################################################################
  
  yhat3 <- data.frame(cbind(yhat2, yhat)) %>% 
    mutate(diff_rep = rephat2-rephat) %>% 
    mutate(diff_dem = demhat2-demhat) %>% 
    mutate(diff_abst = absthat2-absthat) 
  
  
  results[i, 1] <- mean(yhat3$diff_rep)
  results[i, 2] <- mean(yhat3$diff_dem)
  results[i, 3] <- mean(yhat3$diff_abst)
  
  print(i)}, error=function(e){cat("Error",
                                   conditionMessage(e), "\n")}) 
  
}
#18% Missing
############### Simulated vote share
glimpse(results)
sum(is.na(results$V1))/length(results$V1)
18/3000

results <- results %>% 
  rename(rep = V1, dem = V2, abst = V3)

glimpse(results)
results <- results[rowSums(is.na(results)) == 0, ]
glimpse(results)

party <- c("Republican", "Democratic", "Abstention")
simul_data <- as.data.frame(matrix(NA, nrow = 3, ncol = 3))
simul_data <- simul_data %>% 
  rename(party = V1, mean = V2, sd = V3)
glimpse(simul_data)
simul_data[, 1] <- party
simul_data[1, 2] <- mean(results$rep)
simul_data[2, 2] <- mean(results$dem)
simul_data[3, 2] <- mean(results$abst)

simul_data[1, 3] <- sd(results$rep)
simul_data[2, 3] <- sd(results$dem)
simul_data[3, 3] <- sd(results$abst)

#One-tailed t-test
simul_data <- simul_data %>%
  mutate(lb = mean - 1.65*sd) %>% 
  mutate(ub = mean + 1.65*sd)

glimpse(simul_data)

fig_sim <- simul_data %>% 
  mutate(party = factor(party, levels = c("Republican",
                                          "Democratic",
                                          "Abstention"),
                        labels = c("Republican",
                                   "Democratic",
                                   "Abstention"))) %>% 
  ggplot(aes()) +
  geom_hline(yintercept = 0, colour = "gray27", lty = 2, linewidth = 0.75) +
  geom_pointrange(aes(x = party, y = mean, ymin = lb,
                      ymax = ub),
                  lwd = 0.6, position = position_dodge(width = .5), shape = 15
  ) + 
  theme_bw() + #ylim(-0.025, 0.02) +
  labs(y = "Difference in Predicted % of Voters") +#, 
  #   caption = "Effect of +1 Standard Deviation Change in Unemployment Rate") +
  xlab(NULL) +
  theme(axis.text.x = element_text(size = 11#, 
                                   #hjust = 1, 
                                   #angle = 15
  ),
  axis.text.y = element_text(size = 12),
  axis.title.y = element_text(size = 12),
  legend.text = element_text(size = 13),
  legend.title = element_text(size = 13),
  legend.position = "bottom",
  legend.direction = "horizontal",
  title = element_text(size = 14)) #+
#ggtitle("Change in Percentage of Voters after 1-SD Shock in House Price")
fig_sim

glimpse(simul_data)
simul_data <- simul_data 

write.csv(simul_data, here("Data",
                           "AME_Bootstrap_CountyLevel_60Shock.RData"))
